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ABSTRACT 

Numerical simulations of turbulent Rayleigh-Benard convection in an ideal gas, using either 
the anelastic approximation or the fully compressible equations, are compared. Theoretically, 
the anelastic approximation is expected to hold in weakly superadiabatic systems with e = 
A.T/Tr ^ 1, where AT denotes the superadiabatic temperature drop over the convective layer 
and Tr the bottom temperature. Using direct numerical simulations, a systematic comparison 
of anelastic and fully compressible convection is carried out. With decreasing superadiabaticity 
e, the fully compressible results are found to converge linearly to the anelastic solution with 
larger density contrasts generally improving the match. We conclude that in many solar and 
planetary applications, where the superadiabaticity is expected to be vanishingly small, results 
obtained with the anelastic approximation are in fact more accurate than fully compressible 
computations, which typically fail to reach small e for numerical reasons. On the other hand, if 
the astrophysical system studied contains e ^ 0(1) regions, such as the solar photosphere, fully 
compressible simulations have the advantage of capturing the full physics. Interestingly, even in 
weakly superadiabatic regions, like the bulk of the solar convection zone, the errors introduced by 
using artificially large values for e for efficiency reasons remain moderate. If quantitative errors of 
the order of 10% are acceptable in such low e regions, our work suggests that fully compressible 
simulations can indeed be computationally more efficient than their anelastic counterparts. 

Subject headings: convection — Earth — planets and satellites: gaseous planets — Sun: interior — 
Turbulence 


1. INTRODUCTION 

Thermal convection is of primary importance in 
astrophysical objects. It carries the heat flow over 
large regions in stellar and planetary interiors, and 
is one of the major sources of mechanical mixing 
in these objects. Furthermore, some of the most 
striking, large-scale features of stars and planets 
are powered by convective motions, such as in¬ 
trinsic dynamo-generated magnetic fields, plate 
tectonics on Earth and possibly also the zonal 
winds on Jupiter and other giant planets (e.g. 
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The convective regions in stellar and planetary 
objects typically feature a non-negligible density 
stratification, and the flows are often subsonic. In 
this paper, we compare two approaches that are 
commonly used for modeling convection in these 
systems numerically—the fully compressible ap¬ 
proach and the so-called anelastic approximation. 
Our goal is to quantify the accuracy and efficiency 
of both methods in a given situation, guiding mod¬ 
elers in making the right choice for their particular 
problem at hand. 

The fully compressible equations are the most 
fundamental equations governing thermal convec¬ 
tion. They can directly be derived from first prin¬ 
ciples of physics, such as mass, energy, and mo¬ 
mentum conservation, equipped with constitutive 
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relations that characterize the fluid. The resulting 
equations are thus very general and encompass the 
full range of physical behavior, from the temporal 
evolution of the convective motions to the prop¬ 
agation of sound waves. On the one hand, this 
allows to study regions such as the outermost lay¬ 
ers of the Sun, where the Mach number, i.e. the 
ratio of convective velocity to the sound speed, 
becomes 0(1). On the other hand, problems arise 
in low Mach number regions where the flow ve¬ 
locities are much slower than the sound speed, 
which is typically the case in the bulk of deep 
stellar and planetary interiors. Even though the 
convective motions in such regions occur on time 
scales which are many orders of magnitude larger 
than the acoustic time scale, standard numerical 
schemes have to explicitly resolve the sound waves 
for stability reasons. This forces modelers to as¬ 
sume artificially large Mach numbers, which re¬ 
duces the differences between the convective and 
acoustic time scales to numerically tractable val- 


ues (e.g. 

Tobias et aZ.|1998 

Brummell et aL|2002 

Kapyla et al. 

to 

o 

1 —‘ 

o 

|. Errors introduced by this 


procedure occur as an unavoidable side-effect in 
the fully compressible framework. Still, most of 
the numerical resources typically go into captur¬ 
ing acoustic wave propagation phenomena, which 
are generally believed to be irrelevant for the in¬ 


vestigated convection dynamics (but see Bogdan 
et aL|[T9M Meakin and Arnet^|2006 ). 


To circumvent the problems arising from the 
numerical stiffness of the fully compressible equa¬ 
tions, different ’’sound-proof’ models, such as the 


low Mach number approach (e.g. 

Majda and 

SethianI |1985| |Bell et a/.||20041 |Almgren et al. 

20061, th 

e pseudo-incompressible apj 

)roximation 

)roximation 

(Durran 

19891 or the anelastic apj 

(Batchelo 

r 1953 Ogura and Phillips||1962 Gough 

19691 Gilman and Glatzmaier 1981 

Lantz and 

Fan 11999 1 have been developed. Instead of pre- 


scribing artificially large Mach numbers, all these 
approaches take the opposite route by considering 
the small Mach number limit of the fully compress¬ 
ible equations. The same time scale disparities 
which make solving the fully compressible equa¬ 
tions numerically challenging are thus exploited 
to substantially simplify the equations. As a re¬ 
sult, the pressure field adapts instantaneously, 
which effectively filters out the sound waves. This 
comes, however, at the price of loosing the ability 


to study regions where the Mach number is not 
small. 


Among the sound-proof approaches described 
above, the anelastic approximation is the one most 
commonly deployed for modeling stellar and plan- 


etary interiors (e.g. 

Glatzmaier and Roberts|1996 

Miesch et aL||20081 

Brun et al. 20111 Jones et al. 

20111. The anelasi 
predicted to hold f 
in which only slight 
from a hydrostatic 

ic equations are theoretically 
or low Mach number systems 
thermodynamic perturbations 
background state occur (e.g. 
;onvective systems, the back- 
cally assumed to be adiabatic. 
Qs are believed to be satisfied 
s of giant planets and in the 
rvection zone, but break down 
parts which feature relatively 

Gough 19691. In 

ground state is typi 
The above conditio 
in the deep interior 
bulk of the solar cor 
in their outermost 
small sound speeds 

(Ulrich||1970l Bahcall and Ul- 

rich|1988 Christensen-Dalsgaard et a/.|1996| Guil- 

lot et al. 20041. The dynamics of these outer layers 


thus cannot be accounted for within the anelas¬ 
tic framework, and modelers are forced to exclude 
them from the simulation domain. The dynamical 
consequences of neglecting these regions remain 
unclear. 


In summary, both approaches have advantages 
and drawbacks. While the fully compressible ap¬ 
proach is the method of choice for modeling 0(1) 
Mach number flows in near-surface regions of stel¬ 
lar objects, the anelastic approximation seems to 
be beneficial in the deep interiors where the Mach 
numbers are usually very small and where the 
thermodynamic state is close to the adiabat. Un¬ 
fortunately, in many astrophysical applications, it 
remains unclear which approach performs best, 
with anelastic and fully compressible models be¬ 
ing used side by side. The main goal of this study 
is thus twofold: First, we aim to quantify and 
compare the errors inherent in modeling turbulent 
convection in both approaches. Secondly, we seek 
to compare their computational efficiency, thereby 
guiding modelers in minimizing the tradeoff be¬ 
tween accuracy and efficiency for any given situa¬ 
tion. 


Perhaps somewhat surprisingly, comparing re¬ 
sults from the anelastic models currently used in 
astro- and geophysics to standard fully compress¬ 
ible simulations is non-trivial. This is because 
the anelastic models usually parameterize the tur¬ 
bulent, subgrid-scale entropy flux, while similar 


2 





















































































turbulence models are typically not used in fully 
compressible models. The popularity of turbu¬ 
lence modeling in the anelastic framework stems 
from the fact that it allows further simplifications 
of the governing equations, which eases the nu¬ 
merical implementation considerably. Typically, 
molecular heat conduction is neglected and re¬ 
placed by an artificial eddy diffusion model that 


represents turbulent mixing of entropy ( 

Gilman 

and Clatzmaier |I98I[ 

Clatzmaier I984| 

Bragin- 

sky and Roberts|1995 

Lantz and Fan|I999|. This 


turbulent entropy diffusion model, however, is not 
mandatory for the actual anelastic approximation, 
and anelastic equations have been formulated that 
do not rely on parameterizations of the subgrid- 
scale transport (Gough 1969). These equations 


have not found widespread use so far. In order to 
provide direct comparability between the anelastic 
and the fully compressible approach, in this study 
we will restrict ourselves to molecular thermal heat 
conduction in both cases. 

While direct comparisons of anelastic and fully 
compressible gravity wave dynamics in stably 
stratified set-ups have been performed in several 


studies (e.g. 

Davies et al.||2003 

Klein et al.||2010 

Brown et al. 

2012 

), the unstable thermal convec- 


tion case considered in this paper has received 
less attention so far. The work of IBerkoff et o,l.\ 


(2010) focussed on linear magneto convection and 


found good agreement between both approaches 
for the weakly superadiabatic case. Subsequently, 


Lecoanet et al. (2014) studied differences between 


temperature and entropy diffusion, while Calkins 


et al. (2014, 2015) focussed on the influence of 


rotation on the onset of anelastic and fully com¬ 
pressible convection. Their linear study identified 
shortcomings of the anelastic equations for rapidly 
rotating, low Prandtl number fluids, where fast 
density oscillations were found to become non- 
negligible. Calkins et al. (2014) conclude that 


fully non-linear studies tracing the validity range 
of the anelastic approximation are crucial in both 
rotating and non-rotating systems, especially in 
the turbulent regime characterized by a broad¬ 
band frequency spectrum. A first step in this 
direction has been taken by [Meakin and Arnett] 
(2007), who compared non-linear anelastic and 


fully compressible simulations of stellar oxygen 
burning. The differing physical processes included 
in each model, however, precluded a one-to-one 



Fig. 1.— Compressible convection is modeled in 
Rayleigh-Benard geometry, i.e. in a Cartesian box 
that is cooled from above and heated from below. 
Gravity g is pointing downward, antiparallel to 
the z-axis. 


comparability of the anelastic and fully compress¬ 
ible influences. 

In this paper, we present the first systematic 
one-to-one comparisons between fully compress¬ 
ible and anelastic numerical simulations of convec¬ 
tion in the fully nonlinear, turbulent regime. As 
a starting point, we neglect important ingredients 
of stellar convection, such as spherical geometry, 
rotation, compositional inhomogeneities, nuclear 
reactions, magnetic fields, penetration and over¬ 
shooting in stably stratified layers, and the corre¬ 
sponding wave-emission. This allows us to quan¬ 
tify the respective errors, as well as the compu¬ 
tational efficiency encountered in both approaches 
in the simplest setup possible. The influences of 
the above physical processes will be investigated 
in future studies. 

The paper is organized as follows: In section 
we start with defining our idealized model, which 
is followed by discussing the fully compressible 
equations along with the anelastic approximation 
in section A brief overview of the applied nu¬ 
merical methods is given in section while a di¬ 
rect comparison of anelastic and fully compressible 
results and the computational efficiencies of both 
approaches are discussed in section [5] Finally, gen¬ 
eral conclusions are drawn in section [H 

2. MODEL 

Fully compressible and anelastic convection in 
an ideal gas are modeled in a plane fluid layer of 
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depth d confined between rigid, horizontal plates, 
as displayed in figure Gravity g = —gz is con¬ 
stant and pointing downward, antiparallel to the 
unit vector z. The fluid is cooled from above and 
heated from below by maintaining a constant, pre¬ 
scribed temperature difference across the layer. 
The remaining boundary conditions are periodic 
in the horizontal directions and no slip at the bot¬ 
tom and the top boundary. The ideal gas is char¬ 
acterized by constant dynamic viscosity fj, = vp, 
heat conductivity k = Cppn and specific heat ca¬ 
pacities at fixed volume and pressure, c„ and Cp. 
The quantities v and k are the kinematic viscos¬ 
ity and the thermal diffusivity, respectively, which 
consistently vary across the fluid layer inversely 
proportional to the density. 

The governing equations for fully compressible 
convection describing the temporal evolution of 
density p, temperature T, pressure p and veloc¬ 
ity V are 

dtp + V ■ (pv) =0, (1) 


p [dtv + (v • V)v] = -Vp - pgz 


-V(V • v) 


( 2 ) 


way 


c„p[9tr+(v V)T]+p(V-v) 

=Cpp [dtT + (v • V)T] - [dtp + (v • V)p] (5) 

by using equations 0 and 0 (details are given in 
appendix]^. As usual, the thermodynamic quan¬ 
tities are decomposed into a time-independent, 
vertically varying, hydrostatic and adiabatic back¬ 
ground state (index A0and a superadiabatic part 
(index S), which is allowed to vary in time and 
space. 


T{t,x) = Ta{z) +Ts{t,x), ( 6 ) 

p{t,x) = pAiz) + ps{t,x.), (7) 

p{t,y:) = pa{z)+ psit,x). ( 8 ) 

While for subadiabatic or stably stratified fluids a 
conductive background state is the better choice, 
the assumption of approximate adiabaticity (i.e. 
isentropy) is justified in superadiabatic regions, 
where convection turbulently mixes the fluid and 
thus homogenizes entropy. The background pro¬ 
file can be derived from hydrostaticity Vp = —pgz 
(i.e. equation ([^ with v = 0) and the thermody¬ 
namic relation 

pTds = CppdT — 6pdp, (9) 


c„p[atr + (v v)T]+p(v-v) = 

1 n 2 


fcV^T -b 2p e^j - 9(V • v)Sij 


p={cp- Cy)pT, 


(3) 

(4) 


with s denoting specific entropy and the dimen¬ 
sionless thermal expansion coefficient being de¬ 
fined as 6p = —{dinp/dlnT). Note that for an 
ideal gas, 6p = 1, see 0. By further assuming 
adiabaticity (i.e. ds = 0)7 the background state is 
found to be 


with t denoting time and etj = ^ {djVi + diVj) be¬ 
ing the strain rate tensor. Equations (I][3) express 
the conservation of mass, momentum and energy, 
respectively, while equation Q is the ideal gas law. 


3. FULLY COMPRESSIBLE AND AN- 
ELASTIC EQUATIONS 

In the following, the anelastic and fully com¬ 
pressible equations are discussed in detail. 


3.1. Reformulation and Non-dimensionali- 
zation 

We begin with reformulating the left-hand-side 
of equation 0 in the more ” anelastic-friendly” 


TA{z)=Tr(l-^z], (10) 

\ Cp Ir J 

/ XC„/(Cp-C„) 

PA{z) = Pr[l - 7frz] , (11) 

/ g \Cp/{cp-c^) 

Pa{^) — (Cp Cy)prTp I 1 --Z 1 , 

\ Cp 1 TP J 

( 12 ) 

where the index r in Ty and pr refers to reference 
values, here defined as the adiabatic values at the 


^We use the term ’’adiabatic” here for constant entropy 
states. More accurately, the background state may be 
called ” isentropic”, which however appears to be less com¬ 
mon in the literature. 
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bottom boundary. Applying the decomposition of 
the thermodynamic variables (6]|8| to equations 
(T]|4| results in 


dt{pA + Ps) + V • [{pA + /Os)v] = 0, 


(13) 


{PA + Ps) [5tV + (v • V)v] = -V(PA +Ps) 


- {pA + ps)g2 + P 


V^v+ -V(V • v) 


(14) 


Cp{pA + Ps) [dt{TA + Ts) + (v • V){Ta + Ts)] 
- [dt{PA+Ps) + (v • y){pA+Ps)] = 

1 2 


ky\TA + Ts) + 2fi 


-^(V-v)4 


(15) 


peradiabatic pressure ps extracts kinetic energy 
from the vertical flows to drive the horizontal mo¬ 
tions (see e.g. Gough|1969 ). The non-dimensional 
thermodynamic quantitie^ then read 


T(t,x) =r^(z) + ers(t,x), (17) 

p{t,x) = pa{z) + epsit,x), (18) 

p{t,x) = pa{z) + eps(,t,x), (19) 


with e = AT/Tj. and the adiabatic background 
state being 


Ta{z) = (1 - Dz ), 

(20) 

Pa{z) = (1 - , 

(21) 

Pa{z) = {1-Dz)'^/^'^-^K 

(22) 


{PA+Ps) = {cp - Cy){pA + Ps){Ta + Ts). (16) 


Within the anelastic approximation, insignifi¬ 
cant terms in the above equations are neglected. 
To judge which terms are negligible, the magni¬ 
tude of each single term needs to be estimated, 
which is best done after a proper rescaling. If not 
stated otherwise, from now on non-dimensional 
variables will be used. All spatial variables are 
scaled with the domain depths d and velocity is 
non-dimensionalized with a convective free-fall ve¬ 
locity Vf = Correspondingly, time 

is scaled with the free-fall time tf = d/vf = 


yjPrd/ {Apg). In choosing these units, we implic¬ 
itly assume that shorter timescales, as for exam¬ 
ple caused by sound waves, play a minor role. The 
scale for temperature T and adiabatic temperature 
Ta is Tr, i.e. the temperature at the bottom of the 
domain, while the superadiabatic temperature dif¬ 
ference AT, as dictated by the thermal boundary 
conditions, scales the superadiabatic temperature 
Ts- Since temperature and density perturbations 
are usually assumed to be closely correlated (see 
e.g. Clayton 19681, the superadiabatic density 


PS is scaled with the approximate superadiabatic 
density jump Ap = prAT/Tr- Consistently, den¬ 
sity p and adiabatic background density pA are 
scaled with pr, which is the adiabatic density at 
the bottom of the fluid layer. While pressure p 
and adiabatic pressure pA are non-dimensionalized 
with (cp — Cy)prTr as suggested by the ideal gas 
law, the appropriate superadiabatic pressure scale 
Apgd can be inferred from the fact that the su- 


Upon dividing equations ( T3p6 ) by PrVj/d, PrP, 
CpPrVfTr/d, and {cp — Cy)prTr, respectively, we ob¬ 
tain 


edtps + V • [{pA + e/9s)v] = 0, 


(23) 


e{pA + eps) [5tv -b (v • V)v] = 


-V 




-pA + eps - {pA + eps)z 


I Pr 

1 9 1 . J 

V Ra 

V^v+ -V(V-v) 

o 


(24) 


{pA + eps) [sdtTs + (v • V)(Ta -b eTs)] 

1 ' 


- eDdtPs -b (v • V) 

1 


1 - 


7 


PA + f^Dps 


V RaPr 


V^Ta + cTs) 


1 Pr 

1/ 

V Ra 



(25) 


PA + 1 pg = {pA + eps){TA + eTs). (26) 

7 

Due to the non-dimensionalization with charac¬ 
teristic scales, all variables and operators are 0(1) 


^Note that as the temperature at the bottom of the domain, 
which is dictated by the boundary conditions, is used to 
scale the temperature, it follows that T{z = 0) = 1. There¬ 
fore, Ts is generally negative for a superadiabatically strat¬ 
ified system as considered here. 
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and the magnitude of each term in equations (17 


26) can be estimated by its prefactor. The non- 
dimensional control parameters e, Ra, Pr, 7 , and 
D occurring in these coefficients are discussed in 
the following section. 


3.2. Control Parameters and Magnitude 
of the Terms 


Seven control parameters determine the fate of 
the convection system governed by ( 20p2 ) and 
( 23p6 ). The superadiabaticity 


AT _ Ap 

Tr Pr 


(27) 


compares the superadiabatic temperature differ¬ 
ence as dictated by the boundary conditions to a 
typical reference temperature that is—as all other 
reference values—evaluated at the bottom. We 
will show later that e constrains the typical Mach 
number M, which is defined as the ratio of a typ¬ 
ical convective free-fall velocity and the speed of 
sound. The Rayleigh number 




gdh 


(28) 


controls the vigor of convection with large p, d, 
and e enhancing and large diffusivities v and n 
reducing the convective vigor. More formally Ra is 
the ratio of the product of the diffusive timescales 
(P jKr cP IVr to the Square of the free-fall timescale 
tj. The Prandtl number 



(29) 


which for the setup chosen here is constant with 
depth, denotes the ratio of momentum diffusivity 
to the thermal diffusivity and therefore is a ma¬ 
terial property. It effectively controls the impor¬ 
tance of inertia in the system, with Pr <C 1 leading 
to strong and Pr ^ 1 leading to weak inertial ef¬ 
fects. The ratio of the heat capacities defines the 
parameter 



while the Dissipation number 


D = 


gd 


(31) 


can be interpreted in several different ways. Its 
name originates from the fact that it constraints 
how much internal energy can be generated by 
viscous dissipation, i.e. D is a measure for the 
significance of viscous heating with 0 < D < 1 . 
This becomes evident from (27) and (30), which 


allow to rearrange the dissipation number to D = 
This alternative formula- 


1 gdAp _ 1 Epot 


7 PrCuAT 7 AEi,_ 

tion reveals that the dissipation number is propor¬ 
tional to the ratio of the available potential energy 
Epot = gdAp, which drives convection, to the typ¬ 
ical internal energy variations AEmt = PrCyAT. 
As viscous heating results from the dissipation of 
convective kinetic energy (for which Ep^t defines 
the upper limit), viscous heating can only signif¬ 
icantly contribute to internal energy variations if 
Epot is of the same order of magnitude as AEint- 
D can also be interpreted to be the inverse adi¬ 
abatic temperature scale height evaluated at the 
bottom boundary. Finally, the dissipation number 
is directly linked to the density contrast x that 
may serve as an alternative parameter. It is de¬ 
fined as the ratio of the adiabatic density at the 
bottom and at the top, 


pT _ Pa{z = 0 ) 

^ Pa{z = I) 


= (1 - D) 


-l/(7-l) 


(32) 


The total mass of the fluid, as determined by 
the initial conditions, and the aspect ratio of the 
periodic box form the last two control parameters. 


The scaled equations ( 23p6 l, which still repre¬ 
sent the full compressible dynamics, can be fur¬ 
ther simplified by noting that the terms in 
equations ( 24p6 ) balance exactly. In the mo¬ 
mentum equation (241, the terms simply rep¬ 
resent the hydrostatic balance of the reference 
state, i.e. —(1 — l/x)/D\/pA — pa^ = 0. Like¬ 
wise, the first two terms in the energy equa¬ 
tion @ PaVzOzTa - (1 - ^/l)vzdzPA = 0 bal¬ 
ance because of ( 20p2 1 . Note that the conduction 
term l/\/ RaPrV^A drops out here because the 
adiabatic temperature gradient is constant in our 
simple model configuratioi0 Finally, in equation 


^For general depth dependent heat conductivities k and adi¬ 
abatic temperature gradients VT^, this term must be re¬ 
tained. It then effectively acts as a heat source or sink 
and drives or hinders convection with the magnitude being 
estimated by the term’s prefactor l/^/RaPr. For astro- 
physical systems that exhibit large Rayleigh numbers this 
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V • (PAV) = 0, 


(37) 


(26), the e° terms pA = PaTa just represent the 


ideal gas law for the reference state. 

By subtracting the terms from ( 24p6 ) and 
dividing by e, we arrive at 

edtps + V • [{pA + eps)v] = 0, (33) 

{pA + eps) [i9tv + (v • V)v] = -Vps 


1 Pr 

r 9 1 . ,1 

V Ra 

V^v-h -V(V • v) 
0 


PS^+\ -5~ V^v+-V(V-v) , (34) 


{pA + eps) [dtTs + (v • V)rs] - Dp^v^ 
- D [dtp, + (v • V)p,] = 


I Pr 

1, 

V Ra 

v-ij - 


(35) 


D PS _Ts PS PS Ts 

l-^PA Ta Pa ^ PaTa 


(36) 


which describe fully compressible convection as 
perturbations from the adiabatic, hydrostatic 
background state. 


3.3. Anelastic Approximation 


The energy conserving anelastic equations, that 
can also be derived more formally by a rigorous 
amplitude expansion (e.g. Gough|1969 Lantz and 


Fan||1999 1, follow from ( 33|36 ) in the limit e —>■ 0, 


resulting in 


magnitude is typically very small and may be comparable 
or even smaller than the magnitude of the terms repre¬ 
senting the usual convective perturbation. For numerical 
simulations that do not reach realistic parameter values, 
the diffusion of adiabatic background temperature, how¬ 
ever, may be of significance. 


PA [9tv + (v V)v] = -Vps 


- Psz + 


1 Pr 

1 

V Ra 

V^v-f -V(V • v) 
0 


(38) 


PA [dtTs + (v • V)Ts] - Dp,v^ 

- D [dtp, + (v • V)p,] = -^2=V^Ts 


/ Pr 


1 . , 1 

2 

+ 2D\ - 

V Ra 

V-ij 

3 (V • v)(5y 

(39) 

D PS 


PS 

(40) 

1 - yPA 

Ta^ 

PA 


Note that for the setup chosen here, the supera- 
diabaticity parameter e drops out of the non- 
dimensional anelastic equationf^ Furthermore, 
the well-known Boussinesq equations describing 
shallow convection follow in the limit Z) —>■ 0. 


In a very simple manner the above equations 
illustrate the neglected physical processes in the 
anelastic and the Boussinesq approximation: The 
continuity equation (33) reveals that by letting 
e —>■ 0, sound waves are effectively filtered out as 
the time derivative term becomes negligible. Fur¬ 
thermore, unpleasant nonlinearities disappear in 
( Mp6 ). In the Boussinesq limit ZI —0, the en¬ 
ergy equation (35) uncovers that pressure loses its 
role in the energy budget, while viscous heating 
can be neglected as the available potential energy 
is much smaller than internal energy variations 
(c.f. 


section 


3.2). Equation (36) further shows 


that the superadiabatic density is directly propor¬ 
tional to the superadiabatic temperature in the 
Boussinesq limit. Finally, the Mach number 


^ ^ ^ y/ Apgd/pr ^ / eP 

Vs ^Jcj,{cp - Cy)TrlCv y 7 - 1 ’ 


based on the free-fall velocity Vf and the speed 
of sound Vs at the bottom of the domain, can be 
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^For the general case that contains the diffusion of back- 
ground temperature, the e-parameter controls the signifi¬ 
cance of this process and is thus retained. 








































estimated from the input parameters. Obviously, 
it is considered to be small in both, the anelas- 
tic and the Boussinesq approximation. Note that 
the Mach number can serve as an alternative con¬ 
trol parameter that replaces e. In solar and giant 
planets’ interiors, where D and 7 — 1 can typically 
be assumed to be 0(1), the square of the Mach 
number is crudely approximated by the superadi- 
abaticity. 


« e, (42) 

which suggests that the anelastic approximation 
holds for M 1 . 


4. NUMERICAL REALIZATION 


The equations governing fully compressible con¬ 
vection ( 33p6 ) are solved on a collocated grid us¬ 
ing second order finite differences and a third order 
upwind method for the advection terms. A semi- 
implicit time stepping scheme based on a third 
order Adams-Bashforth / backward-difference for¬ 


mula (AB3/BDF3) is applied (e.g. Boyd 2001 
Peyret| 2002[ ). All terms except for the vertical 
diffusion terms are treated explicitly. 

The anelastic simulations that will be presented 
in this paper are performed with an anelastic 
code, which is a modihed version of the Boussi¬ 


nesq code by Stellmach and Hansen (2008). It 


uses a mixed pseudo-spectral fourth order finite- 
difference discretization of the spatial deriva¬ 
tives and an AB3/BDF3 time integration scheme, 
which treats all linear terms implicitly. Instead 
of using ( 33][^ ) directly, for numerical reasons it 
turns out to be benehcial to use an equivalent 
formulation based on entropy rather than temper¬ 
ature. The relevant equations ( C5||C7 ) are derived 
in detail in appendix [Cj 


5. RESULTS 


In this section results from a suite of anelas¬ 
tic and fully compressible direct numerical simu¬ 
lations (DNS) are presented in order to test the 
accuracy and efficiency of both approaches in the 
fully nonlinear regime of convection. 


Equations ( 33p6 ) are solved for various su- 
peradiabaticities (0 < e < 0.4), density contrasts 
(2.72 « exp(l) < X < exp(3) « 20.1 or analo¬ 
gously 0.49 < D < 0.86) and Rayleigh numbers 


e 

Ra 

X 

Resolution 

^run 

Re 

0 

HF” 

2.72 

144^ X 129 

519 

25.0 

0.01 

10 ^ 

2.72 

128^ 

146 

25.0 

0.05 

10 ^ 

2.72 

128^ 

326 

25.6 

0.1 

10‘^ 

2.72 

128^ 

463 

26.3 

0.15 

10 ^ 

2.72 

128^ 

148 

27.0 

0.2 

10 ^ 

2.72 

128^ 

151 

27.7 

0.25 

10 ^ 

2.72 

128^ 

77.8 

28.4 

0.3 

10 ^ 

2.72 

128^ 

185 

29.1 

0.35 

10 ^ 

2.72 

128^ 

92.0 

30.0 

0.4 

10 ^ 

2.72 

128^ 

248 

30.9 

0 

HF” 

2.72 

144^ X 129 

4180 

99.7 

0.1 

10*5 

2.72 

128^ 

3519 

102 

0.2 

10*5 

2.72 

128^ 

2990 

104 

0.3 

10*5 

2.72 

128^ 

2937 

107 

0.4 

10 ^ 

2.72 

128^ 

4260 

111 

0 

10 “ 

2.72 

192^ X 193 

3003 

316 

0.1 

10 “ 

2.72 

1923 

1364 

322 

0.2 

10 “ 

2.72 

1923 

1063 

330 

0.3 

10 “ 

2.72 

1923 

2384 

339 

0.4 

10 “ 

2.72 

1923 

2041 

350 

0 

10 '^ 

2.72 

288^ X 257 

1913 

954 

0.1 

10 " 

2.72 

2563 

1373 

973 

0.2 

10 " 

2.72 

2563 

1262 

991 

0.3 

10 " 

2.72 

2563 

1878 

1016 

0.4 

10 " 

2.72 

2563 

1066 

1050 

0 

10 “ 

4.48 

192^ X 193 

2744 

300 

0.1 

10 “ 

4.48 

1923 

1205 

313 

0 

10 “ 

7.39 

192^ X 193 

2535 

294 

0.1 

10 “ 

7.39 

1923 

1422 

299 

0 

10 “ 

12.18 

192^ X 193 

2811 

280 

0.1 

10 “ 

12.18 

1923 

1347 

286 

0 

10 “ 

20.1 

1922 X 193 

2945 

269 

0.1 

10 “ 

20.1 

1923 

1464 

271 


Table 1: Overview of the simulations carried out 
for this study, with Pr = 0.7 and 7 = 5/3 apply¬ 
ing to all simulations. The horizontal dimensions 
of the simulation domain are lx = ly = “2(1, result¬ 
ing in an aspect ratio of two. While the spatial 
resolution is given in the number of x, y, and z 
grid points, trun denotes the run time measured 
in free-fall times, and Re = Vrmsj\/P t jRa is 
the approximated Reynolds number, with the non- 
dimensional root-mean-square velocity Vrms being 
dehned in equation ( [47| ). While all Ra =10"* cases 
result in stationary solutions, the remaining sim¬ 
ulations stay time-dependent. 

























Fig. 2.— Typical volume renderings of the superadiabatic temperature Tg (a) and vertical velocity Vz (b) for 
an anelastic simulation run that reached statistical equilibrium. Red colors denote warm, buoyant material 
and positive Vz, blue signifies cold fluid and negative Vz, and yellow structures refer to intermediate values 
of Tg and Vz- The corresponding parameters are e = 0, Ra — 10^, Pr — 0.7, 7 = 5/3 and x = exp(l) « 2.72. 
Corresponding snapshots taken from numerical simulations of fully compressible convection look qualitatively 
the same and cannot be visually distinguished from the displayed example. A stereoscopic 3-d version of 
these volume renderings, which reflects the full 3-d structures when wearing red-cyan filter glasses, is shown 
in figure in the appendix. 


(10^ <Ra< 10^). See Table 0 for an overview of 
all simulations. The simulation runs with e = 0 
are carried out with the anelastic code, while the 


is the polytropic index, determines the total mass. 
Note that the polytropic index n is often used as 
an alternative parameter to the superadiabaticity 


remaining simulations are executed with our inde- 

e (e.g. 

Brummell et al. 

1996 

1998 

Berkoff et al. 

pendent code for fully compressible convection as 

20101 . 

For X ~ 2.72 and 0.1 < e < 0.4, which are 


trol parameters are kept constant for all simula¬ 
tions, with the Prandtl nnmber set to Pr = 0.7 
and the ratio of specific heats chosen to repre¬ 
sent a monoatomic ideal gas, 7 = §■ The hor¬ 
izontal dimensions of the simulation domain are 
lx = ly = 2 (i, resulting in an aspect ratio of two. 
The total mass of the fluid is determined by the 
initial state, for which we choose the hydrostatic, 
conductive solution with 


typical parameters for this study, the polytropic 
index varies within the range 1.07 > n > 0.37. 

To give the reader a feeling for the level of tur¬ 
bulence reached in our simulations, figureshows 
a typical snapshot of an anelastic simulation run 
that reached statistical equilibrinm. Correspond¬ 
ing snapshots taken from numerical simulations 
of fully compressible convection look qualitatively 
similar. 


T{t = 0,z) =TA + eTs{t = 0) 

= [l-{D + e)z]. (43) 

The integral over the corresponding initial density 
distribution 

p{t = 0 , z) =pA + eps{t = 0 ) 

= [l-{D + e)z]\ (44) 

where 


n = 


7 D 
7 — 11 ? -I- e 


- 1 


(45) 


5.1. Comparison of fully compressible and 
anelastic results 

Global diagnostic quantities can provide a first 
impression as to what extent the anelastic approx¬ 
imation holds. Initially, we vary e and Ra, while 
keeping x = exp(l) « 2.72 constant. Covering 
several orders of magnitude in Rayleigh number 
Ra, figure shows three different global diagnos¬ 
tics plotted against the superadiabaticity param¬ 
eter e. From left to right, the graphs in the top 
row display the heat flux in terms of the Nusselt 
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Fig. 3.— Global diagnostic quantities are plotted against the superadiabaticity parameter e. From left to 
right, the graphs in the top row display the heat flux in terms of a Nusselt number Nu, the root mean square 
velocity Vrms, and the kinetic energy density Ekin- The bottom row shows the same quantities normalized to 
the corresponding anelastic values (e = 0). For all Rayleigh numbers, the fully compressible results converge 
to the anelastic values for e —>■ 0. For large Rayleigh numbers Ra and superadiabaticities e smaller than 0.3, 
the outputs from compressible convection deviate by no more than 30% from the associated anelastic values. 
The error bars given for the Nusselt numbers are estimates based on the difference between temporally 
averaged Nusselt numbers computed at the top and bottom boundary. In all cases, x = 2.72. 


10 
















number 


Nu = ’ 

(46) 

the root-mean-square velocity 


Vrms ~ J 

(47) 

and the kinetic energy density 


Ekin = 2 {P'^ )t^v ’ 

(48) 


where brackets (...) denote temporal (index t), vol¬ 
ume (index v) and/or horizontal (index h) aver¬ 
ages, while an overbar implies both temporal and 
horizontal averaging, .7. = (...)( ^. The bottom 
row shows the same quantities normalized by the 
corresponding anelastic values 

^kin * 

As expected, the fully compressible cases con¬ 
verge to the anelastic results as e is decreased. 
From theoretical considerations (cf. section]^, we 
expect the convergence to be linear in e, a trend 
which is most clearly seen for Ra = 10^, where the 
flow field is stationary. For the fluctuating solu¬ 
tions encountered at larger Ra, higher order terms 
appear to contribute considerably to the dynam¬ 
ics even for moderate superadiabaticities e > 0.3, 
where the linear scaling is observed to break down. 

Another important question is how the ratio of 
fully compressible and anelastic results scales with 
Rayleigh number for a fixed value of e. Interest¬ 
ingly, while both and Ekin/ElY° de¬ 

crease with Ra, the relative heat flux Nu/Nu'^^^ 
increases, without showing any sign of conver¬ 
gence over the range of Rayleigh numbers stud¬ 
ied. Whether Nu/Nu^^^ converges to a finite 
value beyond Ra = 10^ or continues to increase 
monotonically is left to future investigations. The 
answer is of great importance for astrophysical 
systems, which typically have Rayleigh numbers 
much larger than those considered here. 

In summary, perhaps the most important con¬ 
clusion to be drawn from figure is that in the 
turbulent, high Rayleigh number regime, the dif¬ 
ference between the fully compressible results and 
the corresponding anelastic values remain moder¬ 
ate in all cases studied. For e < 0.3, the relative 
deviations are less than 10%, 20%, and 30% for 
e = O.I, e = 0.2, and e = 0.3, respectively. As a 
crude rule of thumb, e thus provides a reasonable 
estimate of the relative error. 


After discussing global diagnostic quantities 
and their variations with e, a more detailed view 
is provided by vertical profiles obtained from hor¬ 
izontal and temporal averages of the solutions. 
Figure shows temperature and density profiles 
for Ra = 10® and various values of e. The top row, 
from left to right, displays plots of total tempera¬ 
ture T = Ta+cTs, superadiabatic temperature Tg 
and comparisons of T with an approximated total 
temperature Tapprox = Ta + Analogously, 

the total density p = pA + eps, superadiabatic 
density ps and a comparison of the approximated 
density Papprox = PA E with p are presented 
in the bottom panels. Note that the profiles of T 
and p for e = 0 in the left column simply represent 
the adiabatic background state. 

Just like the global diagnostic quantities, the 
profiles for T, Tg, p, and ps obtained from com¬ 
pressible convection simulations converge to those 
of the anelastic case as e decreases. When compar¬ 
ing the approximated temperature profiles Tapprox 
with the corresponding T, no difference can be 
seen by eye for e < 0.2, while for e = 0.4 the 
deviations lie within a few percent. Possible de¬ 
ficiencies of the anelastic approximation become 
evident when comparing p and p approx- While the 
e < 0.2 cases still fit very well, larger superadia¬ 
baticities seem to be problematic especially near 
the top boundary, where the deviations almost 
reach 100% for the e = 0.4 case. 

Further depth profiles are shown for the veloc¬ 
ity field in figure The panels from left to right 
show the time averaged standard deviation of the 
horizontal velocity , the time averaged stan¬ 
dard deviation of the vertical velocity and the 
time averaged skewness of the vertical velocity 7 ^^, 
where the definitions 


ax{z)=(^sJ[X-{X)^fY 


lx{z) 




(49) 

(50) 


are used for the time averaged standard deriva¬ 
tion and skewness of a quantity X. In agreement 
with the results presented above, all profiles ob¬ 
tained from compressible convection simulations 
converge to those of the anelastic case as e de¬ 
creases. 

Finally, we focus on the influence of the density 
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Temperature T Superadiabatic temp. Ts T and Tapprox 



Fig. 4.— Depth profiles of temperatures and densities for various superadiabaticities e for the Ra = 10® 
case. The top row, from left to right, displays plots of total temperatures T = + efs, superadiabatic 

temperatures Tg, and comparisons of T with approximated total temperatures Tapprox = Ta + . These 

are reconstructed from the anelastic limit case by extrapolating to finite superadiabaticities e. Analogously, 
in the bottom row the total densities p = pA + £ps, superadiabatic densities ps, and the approximated 
densities papprox = PA + in comparison with p are presented in the panels from left to right. The 

profiles for T, Tg, p, and ps obtained from compressible convection simulations converge to those of the 
anelastic case as e decreases. No deviations between the approximated temperature profiles Tapprox and the 
corresponding unapproximated ones T can be seen by eye for e < 0.2, while for e = 0.4 the deviations lie 
within a few percent. Deficiencies of the anelastic approximation become evident when comparing p and 
Papprox- While the e < 0.2 cases still fit very well, larger superadiabaticities seem to be problematic especially 
near the top boundary. 
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Fig. 5.— The panels from left to right show the standard deviation of the horizontal velocity the 
standard deviation of the vertical velocity and the skewness of the vertical velocity All prohles 
obtained from compressible convection simulations converge against those of the anelastic case (e = 0) as e 
decreases. 


Dissipation number D Dissipation number D Dissipation number D 



Fig. 6.— Global diagnostic quantities obtained from fully compressible (e = 0.1) and anelastic (e = 0) 
simulations {Ra = 10®) are plotted against the density contrast x (and Dissipation number D) in the top 
row. From left to right, the graphs display the heat flux in terms of a Nusselt number Nu, root mean 
square velocity Vrms, and kinetic energy density The bottom row shows the ratio of the respective 

quantities To guide the eye, a linear fit is also plotted revealing that the relative differences 

of compressible and anelastic outputs generally decrease with increasing density contrasts. The error bars 
given for the ratio of the Nusselt number result from differences in the time-averaged bottom and top Nusselt 
numbers. 
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Fig. 7.— Depth profiles of compressible (e = 0.1) and anelastic (e = 0) convection are shown for different 
density contrasts x with constant Ra = 10®. The panels from left to right show the mean superadiabatic 
temperature Ts, the standard deviation of the superadiabatic temperature arst and the standard deviation 
of the horizontal velocity In the lower half of the fluid container, which covers the largest amount of the 
fluid’s mass, the deviations between the fully compressible and the anelastic case decrease with increasing 
density contrast, while no distinct trend can be seen in the top part. 


contrast x on the accuracy of the anelastic approx¬ 
imation by comparing anelastic (e = 0) and fully 
compressible (e = 0.1) simulations for Ra = 10® 
and 2.72 < x < 20.1. In analogy to figurej^ figure 
[^displays the ratio of compressible and anelastic 
diagnostic outputs. Shown are the Nusselt num¬ 
ber, root-mean-square velocity and kinetic energy 
density for various density contrasts. Interest¬ 
ingly, the differences between anelastic and fully 
compressible diagnostics decrease with increasing 
density contrast. 

The plots shown in figure provide further 
insight into this nonlinear effect. The panel on 
the left shows depth profiles of the superadiabatic 
temperature. They reveal that the magnitude of 
Ts in the lower part of the fluid container, which 
contains the bulk of the fluid’s mass, decreases 
as the density contrast increases. This implies 
that increasing the density contrast x effectively 
decreases the superadiabatic perturbations in the 
bulk region. As the accuracy of the anelastic ap¬ 
proximation is proportional to the relative magni¬ 
tude of the superadiabatic perturbations, increas¬ 
ing the density contrast effectively increases the 
precision of the anelastic equations. 

This view is further supported by the two pan¬ 
els on the right of figure that display the depth 
profiles of the standard deviation of the superadia¬ 
batic temperature ots and the standard deviation 


of the horizontal velocity (jy^. Both exhibit pro¬ 
nounced maxima near the boundaries that mark 
the edges of the thermal and viscous boundary 
layers. The case with the biggest density con¬ 
trast X = 20.1 clearly reveals that the differences 
between anelastic and compressible profiles are 
largest near the location of the top maxima, in 
accordance with the relatively large magnitude of 
the superadiabatic temperature within the upper 
thermal boundary layer. Near the bottom bound¬ 
ary, where the superadiabatic perturbations are 
small, the anelastic curve cannot be visually dis¬ 
tinguished from the compressible profile. 


5.2. Computational Efficiency 


After discussing the quantitative differences be¬ 
tween fully compressible and anelastic results, we 
briefly turn to the question of which approach is 
computationally more efficient. 


The main computational benefit of employing 
the anelastic approximation is that sound waves 
are filtered out and thus do not need to be resolved 


(e.g. |Gough 1969 Glatzmaier 1984 Lantz and Fan 
19991). Assuming that the timestep length in sim¬ 


ulations of anelastic convection is constrained by 
the free-fall velocity u/, while in the fully com¬ 
pressible case the limit is set by the sound speed 
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Vst we expect 


\±compress 

‘-^^max _ ^ 

\4-a7ielastic 9 , 
^^max 


I eP 

7-1 


7-1 


(51) 


where (41) and (32) have been used. The Mach 


number M, or, alternatively, for fixed x snd 7 , the 
superadiabaticity e is therefore expected to control 
the time step ratio. 

In order to check the validity of the above esti¬ 
mate, figure compares the Mach numbers 


M" 


max(^|v| j 

Vs{z = 0 ) 


(52) 


obtained from numerical simulations with the the¬ 
oretical prediction (51). While the maximum is 
taken over the entire computational domain, the 
sound speed is evaluated at the bottom boundary, 
where it becomes largest. Both e and x 9-^6 var¬ 
ied for a mono-atomic ideal gas with 7 = 5/3. As 
expected, the free-fall estimate (51) slightly over¬ 
estimates the real Mach numbers, but overall the 
behavior is captured reasonably well. 


In general, e is the most important parame¬ 
ter controlling the relative efficiency of both ap¬ 
proaches. As figure 1^ shows, for e = 0.1, Mach 
numbers around 0.3 are reached, such that anelas- 
tic codes can use time steps which are roughly 
three times larger than those possible in fully com¬ 
pressible simulations. This however is not guaran¬ 
teed to result in real computational savings, as 
individual time steps tend to be more costly in 
anelastic simulations. In contrast to the fully com¬ 
pressible case, the pressure field adapts instanta¬ 
neously and is governed by an elliptic equation, 
which complicates the time stepping procedure. 
Although efficient solution techniques are read¬ 
ily available from the extensive literature on in¬ 
compressible computational fluid dynamics, it is 
not unreasonable to assume that the additional 
costs slow down the computation of a single time 
step by a factor of three, such that both ap¬ 
proaches reach similar efficiency. Indeed, we ex¬ 
perienced that some of our compressible simula¬ 
tions at e = O.I were computationally more effi¬ 
cient than their anelastic counterparts. 


The above conclusions have been drawn solely 
from simulations of turbulent, compressible Ray- 


Dissipation number D 
0.49 0.74 0.81 0.86 



Fig. 8 .— Plot of the theoretically and numerically 
found Mach numbers M for a monoatomic ideal 
gas, i.e. 7 = 5/3. Numerical M are obtained from 
simulations at Ra = 10® and Pr = 0.7 and show 
good agreement with theoretical predictions that 
are generally slightly larger. While the left panel 
(red) shows a plot of M against superadiabaticity 
e for a constant density contrast of y = 2.72, the 
right panel (blue) displays M against y (or alter¬ 
natively n) for a fixed e = 0.1. The Mach number 
can be interpreted to be the ratio of the respec¬ 
tive maximum timestep lengths of compressible 
and anelastic convection 

which allows to compare the theoretical efficiency 
of anelastic and compressible numerical codes. 


leigh-Benard convection. Effects occurring from 
strong rotation or magnetic fields may alter the 
picture significantly. For example, if thin bound¬ 
ary layers need to be resolved, the time step re¬ 
striction arising from the sound waves may be¬ 
come prohibitive in fully compressible simulations. 
This situation could arise for example in planetary 
cores, where the presence of rigid walls combined 
with the rapid rotation generates very thin Ek- 
man boundary layers, which nevertheless appear 
to be dynamically active (Stellmach et al. 2014) 
and thus need to be accounted for. 

We finish this section by noting that more de¬ 
tailed efficiency comparisons are beyond the scope 
of this paper. Much will depend on the numerical 
algorithms employed, on the degree to which the 
codes are tuned to the machine they run on, and 
on the architecture of the computer itself. 
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6. CONCLUSIONS 


This paper has presented the first one-to-one 
comparison of anelastic and fully compressible tur¬ 
bulent convection. Our goal was to quantify the 
accuracy and efficiency of both methods for the 
simple test case of turbulent Rayleigh-Benard con¬ 
vection in an ideal gas. 


The relation between the anelastic and fully 
compressible equations has been carved out in 
detail, without invoking subgrid-scale turbulence 
modeling at any stage. The anelastic approxi¬ 
mation is expected to hold in the limit of small 
superadiabaticity e, such that the Mach number 
M ~ ^/e remains small. We have shown that the 
fully compressible equations can be manipulated 
into a particular non-dimensional form, which con¬ 
sists of terms representing the anelastic dynamics 
plus 0{e) correction terms that guarantee the fully 
compressible physics. In the limit e —?► 0 these cor¬ 
rection terms vanish and the usual anelastic equa¬ 
tions, as rigorously derived by formal amplitude 
expansions in previous works (e.g. Gough 1969 


Lantz and Fan 1999), are recovered. Our approach 


helps to make the relation between the anelastic 
and fully compressible equations fully transparent, 
and also reveals that the familiar Boussinesq equa¬ 
tions result in the double limit e —)■ 0,1? —>■ 0, 
where D = gd/{cpTj.) is the dissipation number. 
The requirement of small Dissipation number is 
equivalent to a shallow convective system with a 
depth that is much smaller than the typical tem¬ 
perature scale height. 

A key aspect of this work was to quantify the 
differences between fully compressible and anelas¬ 
tic results. Therefore a suite of anelastic and 
fully compressible numerical simulations of ther¬ 
mal convection has been carried out. We com¬ 
pared global diagnostic quantities as well as depth 
profiles of the most important statistical moments 
of thermodynamic variables and velocities. All 
simulations reveal a coherent picture, showing 
that the fully compressible results converge to the 
anelastic ones with decreasing e. The relative de¬ 
viation between both cases was found to be ap¬ 
proximately equal to the superadiabaticity e, indi¬ 
cating linear convergence as predicted by theory. 
For e > 0.3 this linear trend is broken and larger 
deviations are encountered. Besides depending on 
the superadiabaticity, the degree to which both 


approaches give consistent results is controlled by 
the density contrast of the system, i.e. the ratio 
of bottom to the top density. Interestingly, due to 
a nonlinear effect, larger density contrasts reduce 
the quantitative differences between anelastic and 
fully compressible models in our simulations. 


A further aspect of this work dealt with the 
comparison of the numerical efficiency of the 
anelastic and the fully compressible approach. 
In cases with e 1 it is usually argued that 
solving the anelastic equations is computationally 
more efficient than solving their fully compressible 
counterparts, because numerically costly sound 
waves are filtered out. While our results gener¬ 
ally confirm this argument, they also show that 
fully compressible models appear to become more 
efficient than anelastic simulations for e > 0(0.1). 

The implications of our results for the simula¬ 
tions of astrophysical flow phenomena might be 
illustrated by considering the specific example of 
solar convection. Standard solar models like, for 
exam ple. Model S ([Christensen-Dalsgaard et al. 


1996 allow to estimate the superadiabaticity 


e{z) = [dzT — {dzT)yi] jd^T as a function of depth 
within the solar interior. The superadiabaticity is 
predicted to be many orders of magnitude smaller 
than one in the lower 99% of the convection zone. 
It only reaches 0(1) values within the outermost 
one per cent, close to the solar photosphere. This 
suggests that results obtained using the anelas¬ 
tic equations are indeed highly accurate in mod¬ 
els excluding the thin outermost layer where the 
approximation breaks down. The dynamical con¬ 
sequences of neglecting this layer, however, need 
further investigation. 


In contrast, the fully compressible approach in 
principle is capable of capturing the relevant phys¬ 
ical processes throughout the entire convection 
zone. This, however, forces modelers to use un¬ 
realistically large values for e in the bulk of the 
convection zone for numerical reasons. An impor¬ 
tant result of our study is that this procedure in¬ 
troduces moderate errors only. Even for e « O.I, 
where fully compressible codes tend to become 
more efficient than anelastic models, the error in 
global diagnostics such as the overall heat trans¬ 
port or the average kinetic energy was found to be 


®Model S is available online at http://astro.phys.au.dk/ 
~ j cd/solar_iiiodels/ 
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of the order of 10%. The impact on the turbulent 
flow statistics was also shown to remain modest. 
In comparison to other sources of errors, arising 
for example from the inability to reach a realistic 
turbulence level in numerical simulations, a ten 
percent error seems tolerable. 


The above conclusions have been drawn from 
numerical simulations that neglect important in¬ 
gredients of stellar convection, such as spher¬ 
ical geometry, rotation, compositional inhomo¬ 
geneities, nuclear reactions, magnetic fields, pene¬ 
tration and overshooting in stably stratified layers, 
the corresponding wave-emission, and of course 
they did not reach the extreme flow conditions of 
the solar interior. While boundary layers play an 
important role in regulating the convection effi- 


mann and Lohse||2000 

Petschel et aL||2013), their 

dominance is less evidei 
involving much higher 

it in a more realistic model 

dayleigh numbers (|Kraich- 

nan 

1962| |Spiegel| |1971| |Crossmann and Lohse 


2000[) and more realistic boundary conditions 
(Brummell et al. 20021. Applying our results to 


the Sun is therefore somewhat speculative and 
the inclusion of additional physical processes in 
future comparative studies is clearly desirable. In 
particular, an issue that might arise in rapidly ro¬ 
tating systems has recently been pointed out by 


Calkins et al. (2014|, who argued that the anelas- 
tic approximation breaks down in the geo- and 
astrophysically relevant case of rapid rotation and 
low Prandtl number. A study similar to the one 
presented here, but including the effects of rapid 
rotation, is needed to resolve this question and is 
currently underway. 
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APPENDIX 

A. VOLUME WORK TERM AND ENERGY EQUATION 


The volume work term p(V • v) in equation (|^ is based on the divergence of the velocity field, which is 
problematic in the anelastic approach. In the anelastic approximation the velocity field is constrained by 
the anelastic continuity equation V • (pav) = 0, which misses the information of the small superadiabatic 
density changes driving convection. Spiegel and Veronis (1960) deal with a related problem when deriving the 
Boussinesq approximation for shallow convection in an ideal gas. They show that p(V • v) is non-negligible, 
as small superadiabatic density variations become important, although the Boussinesq continuity equation 
requires incompressibility, i.e. V • v = 0. 


Following their procedure, the volume work term in the fully compressible energy equation ([^ can be 
reformulated by using the full continuity equation 0 and the ideal gas law Q. In order to do this, both 
(dimensional) equations are reorganized as follows. 


+ V • (pv) = 0 => V • V = - i [dtp + (v • V)p], 

p={cp- Cy)pT =>dp= - - ' - ^dT. 

{Cp-Cy)T T 


(Al) 

(A2) 


Using the above expressions, the volume work term in 0 can be formulated in terms of temperature and 
pressure, rather than with the divergence of the velocity field. 


piy . v) = - - [dtp + (v • V)p] 

P 

={cp - Cy)p [dtT + (v • V)r] - [dtP + (v • V)p]. 


(A3) 


The anelastic expression for the left hand-side of the energy equation then can be derived by decomposing 
the thermodynamic variables in an adiabatic and a superadiabatic part, as described in section [3.1[ and by 
neglecting all terms involving nonlinearities of variables denoting the superadiabatic part, 

(cp - Cy)p [dtT + (v • V)r] - [dtp + (v • V)p] 

={cp - Cy){pA + ps) [dtTs + (v • V)(rA -b Ts)] - [dtps + (v • V)(pA +Ps)] 
a“=’*‘iC(cp - Cy)pA [dtTs + (v • V)(Ta -b Ts)] + (cp - Cy)d:,TAVzPS - [dtps + (v • V)(PA +Ps)] • (A4) 


The anelastic energy equation then results in 


CpPA [dtTs -b (v • V)Ts] -b CpdzTAVzPs - [dtPs + (v • V)ps] 

n 2 


=kV^{TA+Ts) + 2p 




(A5) 


In contrast, Rogers and Glatzmaier (2005) and Glatzmaier (2014) derive a different version of the volume 
work term that is not equivalent to ours. Instead of using the full continuity equation and the full ideal 
gas law and finally making the anelastic approximation (i.e. neglecting all terms that are proportional 
to nonlinearities of superadiabatic thermodynamic quantities), they apply the anelastic versions of both 
equations. 


V • ( pav ) = 0 


^ dzPA 

V • V =- Vz 

PA 


PA 


PS_ 

PA 


Ts_ 

Ta 


PS 


(cp - Cy){TAPs + PaTs)- 


(A6) 

(A7) 
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which, when utilized for the volume work term, directly reveals their anelastic version 


p(V ■ v) = - {pA + 


PA 


— (Cp c^) 


Ta. 


Ta9zPa H- dzPAPs + QzPaTs 

PA 


(A8) 


The anelastic energy equation from the Rogers and Glatzmaier (20051; Glatzmaier (20141 point of view, then 
results in 


CypA [dtTs + (v • V)Ts] - (cp - Cy)dzPAVzTs =k'V^(TA + Ts) + 2p 


- ^(v •v)4 


1 2 


(A9) 


This anelastic energy equation differs from ours (A31 and essentially involves the superadiabatic temperature 


as the only time-dependent thermodynamic variable. This handy formulation, however, has the disadvan¬ 
tage that corresponding numerical simulations carried out by us neither showed conservation of energy nor 
matched results with fully compressible numerical simulations. 


B. EQUATIONS GOVERNING FULLY COMPRESSIBLE CONVECTION IN ENTROPY 
FORMULATION 

For some applications the entropy formulation of the energy equation (|^, which, at this point, is given 
in dimensional form 


pT -k (v • V)s] = fcV^T -f 2p 




(Bl) 


might be favorable. It can be derived by applying the (dimensional) thermodynamic relation for entropy 

pTds = CppdT — 6pdp. (B2) 


When assuming (5p = 1, as valid for an ideal gas, the integration of equation (B2) reveals the integrated form 


of the thermodynamic relation for entropy, which directly relates entropy to temperature and pressure 

s Sy — Cp In -p— (cp Cy') In , 

1 J’ Pj’ 


(B3) 


where Sy is the reference entropy evaluated at the bottom of the domain. Decomposing the thermodynamic 
variables as done in section |3.1| and exploiting that the dimensional adiabatic background entropy profile 
reads 


- Cp In (Cp Cy') In -(- Sy - Sy 

1J’ P'P 


(B4) 


allows for the reformulation of the whole set of governing equations (T]|4) in terms of entropy instead of 
temperature. The non-dimensional forms of the entropy formulation of the energy equation (iBlI) and the 


thermodynamic relation for entropy (B3) can be derived by using (B4) and As = CpAT/Ty to scale the 
superadiabatic entropy, 


{pA + eps){TA -I- eTs) [9tss -I- (v • V)ss] = 


1 


\/RaPr 


1 Pr 

1, 

V Ra 

Cij — 


1 

Ss = - 


In I 1 -I- e 


Ta 


l-i)ln 1 + ^^ 
7/ V ^-^PA^ 


(B5) 

(B6) 
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C. EQUATIONS GOVERNING ANELASTIC CONVECTION IN ENTROPY FORMULA¬ 
TION 

The anelastic equations in entropy formulation appear from ( 33p4 ), (361 and ( B5|B6 1 in the limit e = 0. 
In order to arrive at the entropy formulation of the anelastic equations the thermodynamic relation for 
entropy (B6) needs to be derived for the limit case e —>■ 0. As lim[ln(l + ax)/x\ = a the integrated form of 

_ £C—>-0^ 

the therrno dynamic relation for entropy in the limit e = 0 results in 

Ts ^PS 
ss = 7^ - D — . 

Ta pa 


(Cl) 


When now replacing Ts in the anelastic (e = 0) ideal gas law (36), the superadiabatic density reads 

D PS 


PS = 


PASS- 


7-17a 

In the following we will express all dynamically varying thermodynamic variables in terms of ss and ps- 

C.l. Pressure and buoyancy term 


(C2) 


Applying the Lantz-Braginsky-Roberts trick (Lantz and Fan 1999 Braginsky and Roberts 1995) to the 
anelastic (e = 0) momentum equation (34) rearranges the pressure and buoyancy terms 

D 


1 y7 PS - 1 „ 

-Vps-z =- Vps - 

PA PA PA 

= -V^- 


Ps 


= -V^ 


PA 

PA 


7-1 paTa 

dzPA , D 
-^PS^ - 


- ssjz 
PS 


Pa 
ssi. 


7-1 paTa 


- Ss 


(C3) 


C.2. Temperature diffusion 


By invoking equation (Cl I, the temperature diffusion term results in 

V2Ts = (tass + D^'^ 


(C4) 


C.3. Governing equations 

By applying equations (C3) and (C4) to ( 33p4 ), (36), (B5) and (Cl) in the limit e = 0 the anelastic 
equations in entropy formulation appear. 


V • (pAv) = 0 

dtv + (v • V)v = —V — 
PA 

PaTa [dtssiy ■ V)ss] = 


Pr 

’"“■"Vito 


y/ RaPr 


V'v+-(V-v) 

s + D^) 

PaJ 


\ / Pr 

1, 


eij - 


(C5) 

(C6) 

(C7) 


The constant temperature boundary conditions can be implemented by using the linearized form of the 
thermodynamic relation for entropy as given by equation (Cl). The entropy formulation of the anelastic 
approximation ( C5||C7 ) can also be derived directly from the anelastic equations in temperature formulation 
(37p0) by using (Cl). Both formulations are fully equivalent. Note that the the temperature diffusion term 


is often neglected and replaced by a parametrized entropy-based large eddy diffusion model (e.g. Glatzmaier 
1984 Lantz and Fan|[T999 Jones et a/.||2011 Gastine et a7|[2M4 ). 
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D. STEREOSCOPIC 3-D ILLUSTRATIONS OF TURBULENT CONVECTION 


Typical volume renderings of the superadiabatic temperature T 5 (a) and vertical velocity (b) for an 
anelastic simulation run that reached statistical equilibrium are shown in figure This stereoscopic 3-d 
version of figure reflects the full 3-d structures when wearing red-cyan filter glasses. 



Fig. 9.— Anaglyph 3-d version of figure that reflects the full 3-d structures when viewing with red-cyan 
filter glasses. 
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